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We present a new method for the calculation of differential distributions directly in Mellin space 
without recourse to the usual momentum-fraction (or z-) space. The method is completely general 
and can be applied to any process. It is based on solving the Integration-by-Parts identities when 
one of the powers of the propagators is an abstract number. The method retains the full dependence 
on the Mellin variable and can be implemented in any program for solving the IBP identities based 
on algebraic elimination, like Laporta. General features of the method are: 1) faster reduction, 2) 
^ I smaller number of master integrals compared to the usual z-space approach and 3) the master inte- 

O ■ grals satisfy difference instead of differential equations. This approach generalizes previous results 

' related to fully inclusive observables like the recently calculated three-loop space-like anomalous 

dimensions and coefficient functions in inclusive DIS to more general processes requiring separate 
, treatment of the various physical cuts. Many possible applications of this method exist, the most 

notable being the direct evaluation of the three-loop time-like splitting functions in QCD. 
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I. INTRODUCTION 



> 

o 

CO I Achieving high precision in theoretical predictions is vital for the success of present and future collider experimental 
programs, as well as for the effective extraction of new physics from experimental data. A significant part of the 
theoretical work related to the experiment requires the evaluation of differential distributions, with most current 
research efforts focusing on the Next-to-Next-to-Leading Order (NNLO) or a higher level of precision. Examples of 
such distributions are the fully inclusivejlllllllllllla and one-particle inclusive l a la DIS, the energy spectrum 
I of hadrons in e+e" collisions Hol pjl fl2l|7 the total partonic cross-section [T^ Il5|r[la| and rapidity distribution 
01 for Higgs and vector boso n lisr production at hadron colliders, Drell-Yan |l9l. l20l. l2ll|. transverse distribution 
Q^l of hadrons at hadron colliders |22l l23l l2^ l25l l2^ or particle spectra in the decays of muon or heavy flavors 
D . |28L I29I I30I I31I Is^ . I33I I . Another important class of distributions that are universal and thus underlay the description 
' of many physical processes includes the space- and time-like splitting functions |3d IstI Ib^ ] , heavy flavor 

matching conditions [39i | and the heavy quark perturbative fragmentation function |40ll4ll l42l|. 

The various distributions can be classified according to the number of kinematical variables they involve. Clearly, 
. , the larger the number of variables, the more complicated the evaluation of a distribution becomes. In this paper we 
■ will restrict our discussion to the case of distributions with a single kinematical variable. This class of distributions 
involves many important examples - some of them still significant open problems - like the three-loop time- like splitting 
functions in QCD. The extension of our discussion to cases with more than one variable will be rather transparent. 

The choice of the most efficient approach to the evaluation of a particular single-scale distribution depends on its 
degree of " inclusiveness" . The fully inclusive observables, like the fully inclusive coefficient functions in DIS 0, allow 
a simplified treatment based on the optical theorem. This is however a rare situation; most distributions of interest 
involve a specific final state, which requires that all contributing physical cuts of the relevant amplitudes be evaluated 
separately. 

The purpose of this paper is to present a conceptually new calculational method of general applicability. As will 
become clear from the subsequent discussion, this method builds a bridge between two very important and seemingly 
unrelated calculational approaches as it provides a new perspective on the calculation of single-scale distributions. 
Moreover, during all stages of calculation this method requires no custom work and utilizes tools, techniques and 
programs that are publicly available and easy to implement in practice. Our method relies heavily on the Integration 
by Parts (IBP) identities It has the important feature of being formulated in terms of variables that are the 

most natural ones for the effective solving of the IBP identities. 

With the above-described applications in mind, let us properly introduce the type of distributions a{z) that we 
will be dealing with in this paper. Such distributions depend on a single kinematical variable z. For example, z can 
be the energy fraction of a parton produced in e~^e~ annihilation. We will assume that this variable is conveniently 
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normalized: < z < 1. The distribution cr is a scalar that is typically of the following form: 

a{z) = [ dPS(™)|M({in} ^ {out})p<5(z - /). (1) 



The factor dPS^™-* in Eq.(^ is related to the phase-space for the m-particle final state; it also contains the measure 
for the virtual integrations (if present). The precise form of this factor depends on the number of particles in the 
initial state. For a single-particle initial state processes with no virtual corrections it reads: 

m 
i=l 

m = j^Y^K<i}~^J)- (2) 

In the case of processes with two particles in the initial state, dPS'™-* has similar structure. It is detailed, for example, 
in [13. 

Typically, expressions like Eq. are UV and infrared divergent and in the following we assume that all divergences 
have been properly regulated by means of dimensional regularization. Besides z, the distribution (j{z) can depend 
on other parameters. Since their presence is irrelevant to our discussion, we will assume in the following that these 
have some fixed values and we will suppress them in our notations. The function / appearing in the argument of the 
(5-function in Eq.(^) is a dimensionless scalar. Its form is specific for each particular observable. 

As a typical example we will consider the evaluation of the single particle inclusive cross-section for massless quark 
production in the decay of a colorless particle V ^ q + X (see also the appendix). Including the corrections up to 
next-to-leading order in the strong coupling and working in terms of bare quantities (i.e. no UV renormalization is 
performed) one has: 



J 1 [dp,] [dp,] |M(") {V ^q + q)\^5(^z 



2pvPq 



Pi 



+ J [dp,] [dp,] [dpg] [M^^HV ^ q + q + g)[H (^z - 



(3) 



In the example above, jM^'^^l^ denotes the terms proportional to a| in the squared matrix element for the process 
V ^ q + X (see Fig.l). Clearly, the first line in Eq.© corresponds to the tree- level (Born) contribution while the 
second and the third lines respectively contain the contributions from the virtual and real-gluon emission corrections 
at order as- On the above example, Eq.QJ stands for any one of the three lines in Eq.Q. 

Perhaps the most elegant approach to date for the evaluation of distributions of the type Eq. ^ was proposed by 
Anastasiou and Melnikov Q and further elaborated upon in Let us recall the salient features of this method. 

One uses the distributional identity: 

2'KiS{x) = 



X + ie X — ie 



to formally replace all (5-functions appearing in Eq.QJ with propagators coinciding with the arguments of the 5- 
functions, i.e. one introduces the invertible mapping P acting only on (5-functions: 



P 



,\{6{x,) 



(4) 



with c an arbitrary function of the propagators. The utility of the mapping Q is that it allows one to treat the 
object on the right-hand side of Eq.Q) with the usual IBP identities |43|. By solving these identities one reduces 
the initial distribution (t(z) to a combination of a small number of irreducible objects. Eventually, one performs 
the inverse mapping P~^, thus expressing u{z) as a linear combination (with simple known coefficients) of a small 
number of well defined master integrals. From the IBP identities it also follows that the master integrals satisfy a 
system of differential equations, which can be solved to obtain their z-dependence. To fully specify the solutions of 
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the differential equations, one has to prescribe the corresponding boundary conditions; these can be extracted from 
an expUcit evaluation of the master integrals in a particular kinematical point like z ~ 1. 

In the phenomenological applications one also needs the Mellin transform ^ of the distribution in question: 

a{n) = [ dzz'^aiz). (5) 

Performing the Mellin transform results in the evaluation of integrals over, typically, combinations of polylogarithms 
and rational functions of z. At present, and certainly for the case of massless distributioris, there exists a ve ry g ood 
understanding of the mapping between the classes of basic functions in z and n spaces 0, E3- In 

the following discussion we will consider the knowledge of cr(z) as equivalent to that of a(n) and vise versa, i.e. we 
will tacitly assume that one can always perform the needed Mellin or inverse Mellin transforms. That is definitely 
true for the massless case. In more complicated situations one may have to resort to numerical methods to perform 
the inverse Mellin transform [soL Isif . In any case, we need not bother about that point here. We will consider our 
problem as solved, as long as we know either a{z) or fj[n). 



II. THE METHOD 



In the present paper we would like to advocate a new approach to the evaluation of the distribution in Eq.lQJ. 
It aims at the direct evaluation of cr(n) without calculating it first in z space as is done at present. Our proposal 
is to explore the obvious possibility that one can integrate over z before performing the phase-space and/or virtual 
integrations: 

a(n) — I dzz"-a{z) 







1 

dzz" / dPS'(™)|M({in} ^ {out}) f6{z - f) 







dF5(")|M({in}^{out})p^. (6) 

We see that as a result of the interchange of the order of integrations the invariant / enters the integrals as 
a propagator raised to power —n. That power, however, should be treated as an abstract parameter that takes 
arbitrary and not fixed integer values. 

By applying the mapping Eq.Q, and interchanging the order of integration as in Eq.(|SJl, one can bring the original 
problem of calculating cr(z) to the following form: 

where V represents the appropriate measure originating from the real and/or virtual integrations. Pi denote the 
propagators originating from the evaluation of the amplitude and Xi are the arguments of all phase-space (5-functions 
(if present). The argument of the (5-function that defines the observed fraction z is denoted by / and the powers a, 
and s are some fixed integers. This way, we have effectively reduced the problem of the calculation of the differential 
cross-section a{z) to the problem of evaluation of functions of the following general form: 

S{ai,...,ap) ^ j pa, ^ pa, - (8) 

The function S appearing in Eq.|(Hl) can also depend on other fixed parameters. 

In principle, scalar terms like the one in Eq.® can be simplified to a minimal set of terms by applying the Integration 
by Parts (IBP) identities Until now, there has been no effective method for solving the IBP identities when one 
(or more) of the powers oi . . . is an abstract parameter (however see 0)- In the following, we present one very 
efficient approach for solving this problem. 



^ Note that usually the Mellin transform is defined through the variable Af = n-|-1,A'^>1. 



4 



III. METHOD FOR SOLVING THE IBP IDENTITIES IN PRESENCE OF AN ABSTRACT POWER 



It is very well known (see for example the book of Smirnov for detailed introduction) that when applied to the 
object S in Eq.||HJ), the IBP identities result in a system of linear homogeneous equations with rational coefficients 
that relate functions S with arguments shifted by ±1 relative to each other. If all a^'s were fixed integers, then by 
successively relating terms that differ with ±1 one can eventually express the original function S{ai, . . . , Op) through 
a linear combination of several, say to, master integrals S'i({ii}), . . . , S'm({*m})- The masters Sj are special cases of 
S{ai, . . . , Gp) with their arguments {HA) taking special values. At present, the most popular method for solving the 
IBP identities is the one of Laporta |5J| . It is based on solving the systems of linear homogeneous equations directly, 
through Gauss elimination. 

Clearly, if one of the parameters ai is not an integer this procedure cannot work, since: 1) with only integer steps 
one cannot relate the initial non-integer parameter to an element S{ai,...,ap) with only integer ai's and 2) the 
number of steps in the Gauss elimination cannot even be specified when one of the parameters is an abstract number. 

In the following, we present a solution to this problem. In order to facilitate our discussion we shall assume that 
ai, . . . , ap_i are integers having some specific values, while the last argument, Op, is an abstract parameter. 

From Eq.(^ it is clear that (T(n) is a sum of a number of terms of the type in Eq.® that have different values of 
their indexes (ai, . . . , Op). It is very important to observe, however, that the difference between any two values that 
the index Op can take is always an integer, i.e. Op — G N. 

This is a crucial observation, which one can use to modify the strategy for solving the IBP identities in the following 
way. First, one relaxes the requirement that the masters must have integer- valued indexes. Second, as we will explain 
in a moment, one can choose all masters in such a way that they all have the same value, say = r ^ N, of their 
last index i.e. the masters are all of the form: 



S{ji, ■ ■ • , jp-i,r) = V 



1 1 
^1 ■ ■ ■ ^p-i •> 

with the same r, and the jVs being fixed integers specific to each master. 

It is indeed possible to arrange that all masters have the same value of the non- integer- valued index Op. That 
follows from the arbitrariness of the value of this parameter (we only assume that it is non- negative) . Since there is 
no preferred value for that index, the IBP system has a sort of translational invariance along the index Cp. One can 
understand this by saying that r and r + fc, where is a fixed integer, are equally arbitrary. Therefore, we can take as 
a reference value for the index Cp the number r which we will consider abstract but having fixed value. Having done 
that, the "translational" invariance along the values of Op is now "broken". Clearly, the value r now plays the role of 
a zero reference point much like the value Cp = in the usual case when all indexes take integer values. Therefore all 
one need to do is to measure in integer units how much the value of the last index of an element S is displaced from 
the reference point r. 

Next, we give a practical recipe of how to implement the above idea. Let us work with the functions ^ B: 



B{ai, . . . , ap_i, /c) ^ J V — 



(9) 



where as the "reference" point for the last index we take the Mellin variable n. 
As follows from Eqns.© and O the distribution cr(n) takes the following form: 

= X! Cai,...^ap_i-B(ai,...,ap_i,fc), (10) 

ai ,. . . ,ap_i ,/c 

where Cai,...,a _i are some known coefficients. To construct the needed algebraic reductions, one first applies the IBP 
identities on a generic monomial of the form: 

(11) 



where all powers are treated as arbitrary parameters. Next, one identifies each term of the form (|ll|l appearing in 
the IBP equations, with the function B{vi, . . . , fp_i, Vp + n), followed by the substitution Up ^ Vp — n. After this 



^ These are essentially the same as the functions S introduced above; the difference is in the notation used for the last argument. 
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manipulation the Mellin variable n is explicitly present as a parameter in the resulting equations. They can be solved 
in any approach available, including the one of Laporta. 

A word of caution: one has to keep in mind that, as follows from Eq.©, the functions B implicitly depend on n. 
Therefore, one should not confuse the integer value k in the last argument of the function B with the absolute power 
of the corresponding propagator 1//, but should think of it as the "distance" - in integer units - from the reference 
power n. 

Finally, one can map all integrals appearing in Ea. HlUI) to the masters obtained from the solving of the just-described 
reduction. This mapping is done in the standard way. 

Next, we explain how one can extract the n dependence of the master integrals. Assume that an element 
B{bi,. .. ,5p_i,0) is a master integral (with 61, . . . ,6p_i some fixed integers). One can inspect the already solved 
IBP reduction and read off from there the result for the element B{bi, . . . , 6p_i, —1). Note that this element differs 
from the master B{bi, . . . , 6p-i, 0) only by the value of the last index. If the element B{bi, . . . , fcp-i, —1) is not a 
masters itself, then it must be a linear combination of the master integrals: 

. . . , 6p_i, -1) = c(n)B(6i, . . . , 6p_i, 0) + G{n). (12) 

Here c(n) is a known, typically not very complicated function, and the term G{n) is a homogeneous linear combination 
of all master integrals, except for the master B{bi, . . . , bp-i, 0). Ea. (|12|) is a first order non-homogeneous difference 
equation of the type F{n + 1) = c{n)F{n) + G{n) (recall Eq.@) for the master B{bi, . . . , 6p_i, 0). Clearly, repeating 
this procedure for each one of the master integrals found in the reduction run, one can derive a complete system 
of difference equations for all the masters. Typically, one observes certain hierarchy among the master integrals; 
the simplest ones satisfy homogeneous equations (i.e. G(n) = 0) that can be solved in terms of F-functions. These 
integrals then comprise the non-homogeneous terms for the equations of other masters, and so on. 

In case the element i3(6i, . . . , 6p_i, — 1) is also a master integral, one should read off from the reduction the 
result for the yet higher term i3(6i, . . . , 6p_i, — 2). One should continue doing this until one reaches an element 
B{bi, . . . , 6p_i, —k) which is not a master itself but all elements B{bi, . . . , 6p_i, — s) with < s < A: are masters. The 
result from the reduction for the element B{bi, . . . , 6p-i, —k) represents a fc-th order difference equation for the master 
integral i3(6i,...,6p_i,0). 

To solve the resulting difference equations one can make use of existing techniques. Such equations were analyzed 
and successfully solved in the course of the evaluation of the three-loop anomalous dimensions in QCD [s^, HI] and 
of the two '3| and three- loop ^ coefficient functions in DIS. In most cases of physical interest the resulting difference 
equations can be solved after expansion in e in terms of harmonic sums or their generalizations. In simpler cases, one 
can even solve these equations in closed form in terms of hypergeometric and/or F-functions. Difference equations as 
a way of calculating master integrals were also used by Laporta [s^l ■ 

Upon solving the system of difference equations for the master integrals, one has achieved a complete extraction 
of the dependence of the masters on the Mellin variable n. The only remaining thing to do is to specify the initial 
conditions for the solutions of the difference equations. Typically, that would be the value of the masters for n = 0. 
This is an important fact. It implies that to completely specify the master, one need to only evaluate integrals that 
are fully integrated over the available phase-space. These are pure numbers that do not depend on the kinematical 
variable n. On the conceptual level, this is placing the evaluation of certain not-completely inclusive observables one 
step closer to the very familiar fully-inclusive case where, thanks to the optical theorem, one can significantly simplify 
the calculations by not considering separately all possible physical cuts. 

Often, one can reduce the number of fixed-n integrals that have to be evaluated by hand. This follows from the 
property of the fixed-n reduction that not all integrals corresponding to initial conditions for the n-dependent masters 
are actually independent. To explore this fact one has to perform a separate fixed-n IBP reduction where n = is taken 
from the very beginning. We have observed in simple one- and two-loop reductions as well as in rather complicated 
three-loop cases that this procedure indeed generates additional relations between the integrals corresponding to the 
initial conditions of the master integrals. One might wonder about the cost of such an additional run. That, however, 
should be of no concern since the fixed-n reduction is much simpler and faster than the general-n run one has to be 
able to perform anyway. The computer load pays off with the elimination of many of the integrals that otherwise have 
to be computed by hand. A fixed-n reduction can also be used as a cross-check of the general-n calculation. This is 
similar to the use of Mincer in the three-loop DIS calculations in [s^ [s^, 13 • 

IV. PARTIAL FRACTIONING 

Consider a case where during the evaluation of the amplitudes one gets a propagator that is not linearly independent 
from the constraint / defining z (see Eq.j?!)). Clearly, that can happen in many ways and such linear dependence 
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might even involve a group of several propagators. However, to simplify our point as much as possible, we will only 
consider a simple situation. Consider the function: 

^,„, = /p.,..,._l_i, a„ 

where . . . stay for powers of other possible propagators. If we were to evaluate this integral in 0-space we would first 
replace / everywhere with z, as is implied by the factor d{z — /). Then the factor 1/(1 — /) becomes just the number 
1/(1 — z) and drops out of the integral. However, when we work in Mellin space the constraint f = z cannot be 
used anymore. If n were some fixed integer, we could have applied partial fractioning n-times and split the linearly 
dependent propagators 1// and 1/(1 — /). For symbolic n, however, that cannot be done and one should again resort 
to solving difference equations. This can be done in the following way. The identity: 

1 111 1 

(1 - /) " (W)7^^7^' 

immediately translates into a difference equation for the function appearing in Ea. H13|l : F{n + 1) = F{n) + G{n). 
This is a simple difference equation with non-homogeneous part given by: 

G(n) = Jvx{...)xj^. (14) 

The integral (|14|l does not contain linearly dependent propagators and can be evaluated by using the procedures 
described previously. 

Another way of eliminating the linear dependence among the propagators is to expand the propagator 1/(1 — /) in 
geometric series. That would completely eliminate this propagator and one would end up with a standard problem 
where the index n is replaced hy n + k, k > (here one can apply the usual reduction since the index n + fc is as 
arbitrary as the index n is). Finally, one would have to sum up the resulting expression over the index fc. 

Our experience shows that the combination of the above methods is sufficient to eliminate the appearance of linearly 
dependent propagators in any situation. 

V. CONCLUDING REMARKS 

The method presented in this paper represents a conceptually new approach for the evaluation of differential 
distributions. At the same time it has the advantage that it shares features with other available methods that 
previously were applied with impressive success. We would like to critically compare our approach with these methods 
and clarify its distinct applicability. 

Our method has the following advantages over the direct momentum space approach First, the solving of the 
IBP reductions is much more efficient and fast. In practical terms that may not be an issue in simple one- or two- loop 
cases, but we have checked that it brings enormous improvement when applied, for example, at three loops. Second, 
the number of master integrals being produced in the course of solving the IBP's is smaller and, third, our method 
only involves computation of integrals that are pure numbers. This is a significantly easier task compared to the 
calculation of the z-dependent initial conditions at fixed values of z. 

We also expect our method to bring new insights into the calculations of distributions with more than one scale. 

Our approach has a few common features with the method used in the inclusive DIS calculations at two- and 
three-loops H i, HI HUll HI • Still our method has wider, in fact completely general, applicability and it relies 
for the formulation and solving of the IBP identities only on well established, multipurpose and publicly available 
methods and software. An important shared feature is the fact that the master integrals in both methods satisfy 
difference equations in terms of the Mellin variable n. In this respect our method benefits greatly from these existing 
developments, since many of the technical tools needed for its practical realization have been already established. 
That includes the methods for solving the difference equations, the current good understanding of the appropriate 
functional bases in n and z spaces and their tabulation. Finally, well established multipurpose software like FORM 
[ssf. Summer fs^ and XSummer '6^ exist and are publicly available. They are capable of effectively dealing with the 
necessary algebraical manipulations. In fact, our method opens up new venues for the application of the techniques 
from the two- and three- loop inclusive DIS. 

As a first application of our method we have rederived the NLO coefficient functions in e+e~ 10] (in fact to all 
orders in e). 

We also have experience with reductions at three-loops, where the advantage of our method becomes apparent. Our 
method ensures better performance of the IBP solving software and produces smaller number of master integrals. 
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APPENDIX A: SIMPLE EXAMPLE 



We present one example which is simple yet it demonstrates all non-trivial features of the method discussed above. 
We consider the evaluation of the coefficient function in the decay of a colorless object to a quark-antiquark pair 
V ^ q + X at order 0:5. One constructs |Mp in the usual way; the relevant diagrams are shown on Fig.l. At leading 
order, the differential observable of interest is: 

aiz)^^^=Sil-z) + 0{as). (Al) 

Note that we have chosen normalization where the coefficient of the (5-function in the leading term is exactly one to 
all orders in e, d — A — 2e. In Mellin space, this corresponds to cr(n) = 1. The virtual corrections at order as produce 
the same type of contributions. 

Through order as, all non-trivial n-dependence of the distribution Ea. ljAl|l originates from the real gluon radiation 
diagrams. To be specific, we are interested in observing the final state massless quark in the reaction V{p) —> 
q{p ^ Q ^ If') + ^ : where the unobserved massless antiquark and gluon carry momenta q and k respectively. We take 

= 1. From the independent momenta p, q and k one can construct five linearly independent scalars that are needed 
to build the IBP reduction. As such we choose: Pi = (p-q)'^, P2 = g^, P3 = fc^, P4 = ip-q-kY, P5 = 2-2p.q-2p.k. 

After we perform the mapping I^J we have to deal with monomials of the type: 



11111 

pai pa2 paa pa4 pas 
^1 ^2 ^3 4 ^5 



(A2) 



Performing the IBP reductions in the Laporta's method implemented in the program AIR l57l | , we find that there 
is a single master integral: 

B{0, 1, 1, 1,0)= / d''qid\26{Pi)5{P2)m)^ (A3) 

By inspecting the results for the element i?(0, 1, 1, 1, —1) from the solved IBP reduction, we derive the following 
difference equation for the only master integral: 

2 + n ~ 2e 

i?(0,l,l,l,-l) = — 5(0,1,1,1,0). 

3 + rt — 3e 

It is trivial to solve this recurrence relation (we modify the notation in an obvious way to make completely transparent 
the n-dependence): 

The initial condition i?(0, 1, 1, 1, 0)(n ~ 0) is defined through Ea. (|A3p after setting n = there. It is a trivial to 
compute number. 

Our work is not quite done with the solving of the IBP reductions and the evaluation of the master integral since 
there also appears the propagator Padd — {q + k)"^ which does not belong to the set Pi, . . . ,P5. This propagator 
results from the right diagram on Fig.l; it is not linearly independent from the set Pi, . . . ,P5 but it can appears 
downstairs together with these propagators. To resolve the situation one has to resort to the partial fractioning 
technique discussed in section Hvl 

Exploiting the constraints implied by the three (5-functions (with arguments P2.3,i), one can easily establish that 
Padd = 1 ~ ^5- In this case one can apply the geometric series trick discussed in the previous section to all terms 
where Padd is present downstairs: 

1 1 °° 

pn pn \ ^ ps 



Padd 1 - ^5 
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FIG. 1: Real emission diagrams contributing to the decay of a colorless particle V{p) q(p — g — fc) + X at 0{as). 

Next one takes the summation outside the integrals; the resulting integrand is of the type in Ea. ljA2|l . With the help 
of the IBP reduction this integral can be reduced to the master integral discussed above. Inserting the explicit form of 
the master, one obtains very simple expression containing only F-functions. The summation over s of this product of 
F-functions can be easily performed and it results again in a product of the same type of functions. Thus, the result 
from the real emission radiation at order as can be easily evaluated in closed form to all orders in e. To use this 
expression in practical applications one has to decompose it in series in e. This expansion can be easily automated 
with the help of the programs Summer |59| and XSummer |60l | . 
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